#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define NUM_NEURONS	3
#define NUM_SYNAPSES	(NUM_NEURONS - 1)
#define NUM_SAMPLES	2

#define WEIGHT_STEP	0.1
#define ACCEPTABLE_ERROR 0.01

#define INITIAL_TEMP	200.0
#define FINAL_TEMP	1.0
#define TEMP_FACTOR	0.99

#define NUM_ANNEALINGS_PER_SAMPLE	100
#define	NUM_SAMPLES_FOR_AVERAGE		100
#define NUM_ITERATIONS_FOR_CONVERGENCE 	(200 * NUM_NEURONS)

double n[NUM_NEURONS];
double b[NUM_NEURONS];
double w[NUM_SYNAPSES];

double positive_synapses_statistics[NUM_SYNAPSES];
double negative_synapses_statistics[NUM_SYNAPSES];
double positive_bias_statistics[NUM_NEURONS];
double negative_bias_statistics[NUM_NEURONS];

double data[NUM_SAMPLES][2];


void
read_data()
{
	data[0][0] = 1.0; data[0][1] =  0.0;
	data[1][0] = 0.0; data[1][1] =  1.0;
}


void
print_neurons()
{
	int i;
	
	for (i = 0; i < NUM_NEURONS; i++)
		printf("%0.2lf, ", n[i]);
		
	printf("\n");
}


void
print_biases()
{
	int i;
	
	for (i = 0; i < NUM_NEURONS; i++)
		printf("%0.2lf, ", b[i]);
		
	printf("\n");
}


void
print_synapses()
{
	int i;
	
	for (i = 0; i < NUM_SYNAPSES; i++)
		printf("%0.2lf, ", w[i]);
		
	printf("\n");
}


void
randomize_neurons()
{
	int i;
	
	for (i = 0; i < NUM_NEURONS; i++)
		n[i] = (double) (rand() % 2);
}


void
randomize_bias()
{
	int i;
	
	for (i = 0; i < NUM_NEURONS; i++)
//		b[i] = ((double) ((double) rand() / (double) RAND_MAX) * 2.0 - 1.0);
		b[i] = 0.0;
}


void
randomize_synapses()
{
	int i;

	for (i = 0; i < NUM_SYNAPSES; i++)
//		w[i] = ((double) ((double) rand() / (double) RAND_MAX) * 2.0 - 1.0);
		w[i] = 0.0;
}


void
build_network()
{
	randomize_neurons();
	randomize_bias();
	randomize_synapses();
}


double
logit_sample(double logit)
{
	double sample;
	
	sample = (double) rand() / (double) RAND_MAX;
	if (sample > logit)
		return 0.0;
	else
		return 1.0;
}


double
logit(double deltaE, double T)
{
	return (1.0 / (1.0 + exp(-(deltaE / T))));
}


double
logit_function(int i, double T)
{
	double deltaE = 0.0;
	
	if (i != 0)
		deltaE += n[i - 1] * w[i - 1];

	if (i != (NUM_NEURONS - 1))
		deltaE += n[i + 1] * w[i];
		
	deltaE += b[i];
	
	return (logit(deltaE, T));
}


void
permute(int *vector, int size)
{
	int i, j, aux;
	
	for (i = size - 1; i >= 1; i--)
	{
		j = rand() % (i + 1);
		aux = vector[i];
		vector[i] = vector[j];
		vector[j] = aux;
	}
}


void
run_hidden_net(double T)
{
	int i;
	int neuron[NUM_NEURONS - 2];
	double logit;
	
	for (i = 0; i < (NUM_NEURONS - 2); i++)
		neuron[i] = i;
	permute(neuron, NUM_NEURONS - 2);
	
	for (i = 1; i < NUM_NEURONS - 1; i++)
	{
		logit = logit_function(neuron[i - 1] + 1, T);
		n[neuron[i - 1] + 1] = logit_sample(logit);
	}
}


void
run_full_net(double T)
{
	int i;
	int neuron[NUM_NEURONS];
	double logit;
	
	for (i = 0; i < NUM_NEURONS; i++)
		neuron[i] = i;
	permute(neuron, NUM_NEURONS);
	
	for (i = 0; i < NUM_NEURONS; i++)
	{
		logit = logit_function(neuron[i], T);
		n[neuron[i]] = logit_sample(logit);
	}
}


void
update_synapses_statistics(double *statistics)
{
	int i;
	
	for (i = 0; i < NUM_SYNAPSES; i++)
		statistics[i] += n[i] * n[i + 1];
}


void
update_neurons_statistics(double *statistics)
{
	int i;
	
	for (i = 0; i < NUM_NEURONS; i++)
		statistics[i] += n[i];
}


void
init_statistics(double *synapses_statistics, double *bias_statistics)
{
	int i;
	
	for (i = 0; i < NUM_SYNAPSES; i++)
		synapses_statistics[i] = 0.0;
	for (i = 0; i < NUM_NEURONS; i++)
		bias_statistics[i] = 0.0;
}


void
sumarize_statistics(double *synapses_statistics, double *bias_statistics)
{
	int i;
	
	for (i = 0; i < NUM_SYNAPSES; i++)
		synapses_statistics[i] /= ((double) NUM_SAMPLES_FOR_AVERAGE * (double) NUM_SAMPLES * (double) NUM_ANNEALINGS_PER_SAMPLE);

	for (i = 0; i < NUM_NEURONS; i++)
		bias_statistics[i] /= ((double) NUM_SAMPLES_FOR_AVERAGE * (double) NUM_SAMPLES * (double) NUM_ANNEALINGS_PER_SAMPLE);
}


double
get_noise_value(double visible_unit_input)
{
	double noise_value;
	
/*	if (visible_unit_input == 0.0)
	{
		if ((rand() % 100) < 5)
			noise_value = 1.0;
		else
			noise_value = 0.0;
	}
	else
	{
		if ((rand() % 100) < 15)
			noise_value = 0.0;
		else
			noise_value = 1.0;
	}
*/
	noise_value = visible_unit_input;
	return (noise_value);
}


void
set_visible_units(int input_sample)
{
	n[0] = get_noise_value(data[input_sample][0]);
	n[NUM_NEURONS - 1] = get_noise_value(data[input_sample][1]);
}


void
positive_phase()
{
	double T;
	int i, j, input_sample;
	
	init_statistics(positive_synapses_statistics, positive_bias_statistics);
	
	for (j = 0; j < NUM_ANNEALINGS_PER_SAMPLE; j++)
	{
		for (input_sample = 0; input_sample < NUM_SAMPLES; input_sample++)
		{
			randomize_neurons();
			set_visible_units(input_sample);
			
			for (T = INITIAL_TEMP; T >= FINAL_TEMP; T = T * TEMP_FACTOR)
				run_hidden_net(T);

			for (i = 0; i < NUM_SAMPLES_FOR_AVERAGE; i++)
			{
				run_hidden_net(FINAL_TEMP);
				update_synapses_statistics(positive_synapses_statistics);
				update_neurons_statistics(positive_bias_statistics);
			}
		}
	}

	sumarize_statistics(positive_synapses_statistics, positive_bias_statistics);
}


void
negative_phase()
{
	double T;
	int i, j;
	
	init_statistics(negative_synapses_statistics, negative_bias_statistics);
		
	for (j = 0; j < NUM_ANNEALINGS_PER_SAMPLE * NUM_SAMPLES; j++)
	{
		randomize_neurons();
		for (T = INITIAL_TEMP; T >= FINAL_TEMP; T = T * TEMP_FACTOR)
			run_full_net(T);
		
		for (i = 0; i < NUM_SAMPLES_FOR_AVERAGE; i++)
		{
			run_full_net(FINAL_TEMP);
			update_synapses_statistics(negative_synapses_statistics);
			update_neurons_statistics(negative_bias_statistics);
		}
	}

	sumarize_statistics(negative_synapses_statistics, negative_bias_statistics);
}


void
synapse_update()
{
	int i;
	
	for (i = 0; i < NUM_SYNAPSES; i++)
	{
		if (fabs(positive_synapses_statistics[i] - negative_synapses_statistics[i]) != 0.0)
			w[i] = w[i] + WEIGHT_STEP * ((positive_synapses_statistics[i] - negative_synapses_statistics[i])) / fabs((positive_synapses_statistics[i] - negative_synapses_statistics[i]));
		fprintf(stderr, "%lf ", w[i]);
	}
}


void
bias_update()
{
	int i;
	
	for (i = 0; i < NUM_NEURONS; i++)
	{
		if (fabs(positive_bias_statistics[i] - negative_bias_statistics[i]) != 0.0)
			b[i] = b[i] + WEIGHT_STEP * ((positive_bias_statistics[i] - negative_bias_statistics[i])) / fabs((positive_bias_statistics[i] - negative_bias_statistics[i]));
		fprintf(stderr, "%lf ", b[i]);
	}
	fprintf(stderr, "\n");
}


int
convergence_test()
{
	int i;
	double error = 0.0;
	
	for (i = 0; i < NUM_SYNAPSES; i++)
		error += (positive_synapses_statistics[i] - negative_synapses_statistics[i]) * (positive_synapses_statistics[i] - negative_synapses_statistics[i]);
		
	for (i = 0; i < NUM_SYNAPSES; i++)
		error += (positive_bias_statistics[i] - negative_bias_statistics[i]) * (positive_bias_statistics[i] - negative_bias_statistics[i]);
		
	error = sqrt(error);
	//fprintf(stderr, "%lf\n", error);
	
	if (error < ACCEPTABLE_ERROR)
		return (1);
	else
		return (0);
}


void
train()
{
	int converged = 0;
	int num_iteractions = 0;
	
	while (/*!converged && */ (num_iteractions < NUM_ITERATIONS_FOR_CONVERGENCE))
	{
		positive_phase();
		negative_phase();
		synapse_update();
		bias_update();
		converged = convergence_test();
		num_iteractions++;
	}
	printf("converged = %d\n\n", converged);
}


void
run()
{
	int i;
	double T;
	
	for (i = 0; i < 20; i++)
	{
		randomize_neurons();
		for (T = INITIAL_TEMP; T >= FINAL_TEMP; T = T * TEMP_FACTOR)
			run_full_net(T);
	
		run_full_net(FINAL_TEMP);
		print_neurons();		
	}
}


int
main()
{
	srand(time(NULL));
	
	build_network();
	read_data();
	train();
	run();
	printf("\n");
	print_synapses();
	printf("\n");
	print_biases();
	
	return (0);
}
